In cheungngo/diffeq: Studies related to differential equations

library(deSolve)
library(Deriv)
library(Ryacas0)

A simple Michaelian System

yini = c(A = 100, As = 0)
kp = 1; km = 0.3

mich = function(t, y, parms) {
  with(as.list(y), {
    dA = km * As - kp * (100-As)
    dAs = kp * (100-As) - km * As
    list(c(dA, dAs))
  })
}

times = seq(0,10,0.05)

out <- ode(y = yini, times = times, func = mich, parms = NULL)

plot(out)

out

deriv_data = function(y,t) {
  dy = y[-1]-y[-length(y)]
  dt = t[2]-t[1]
  dydt = dy/dt
  return(dydt)
}

# Solving the equation graphically

dydt_As = deriv_data(out[,3],out[,1])
dydt_A = deriv_data(out[,2],out[,1])
plot(dydt_As~out[-1,3], xlim = c(0,100))
lines(dydt_A~out[-1,3])

Finding the solution graphically

FR = function(As) {kp * (100-As)}
BR = function(As) {km * As}
Solve(kp * (100-As) == km * As, As)
FR(100/1.3)
BR(100/1.3)

plot((1:100)/100,FR(1:100));lines((1:100)/100, BR(1:100))

What if Kp is a function of stimulus [S]

yini = c(A = 100, As = 0)
km = 0.5

mich = function(t, y, parms) {
  with(as.list(y), {
    dA = km * As - kp * (100-As)
    dAs = kp * (100-As) - km * As
    list(c(dA, dAs))
  })
}

times = seq(0,10,0.05)

As_steady = c()
for (i in 1:100) {
  kp = i
  out <- ode(y = yini, times = times, func = mich, parms = NULL)
  As_steady = c(As_steady,out[200,3])
}

plot(1:100,As_steady)

Michaelian system with linear feadback kf * A

yini = c(A = 100, As = 0)
kp = 1; km = 0.3; kf = 0.5


mich = function(t, y, parms) {
  with(as.list(y), {
    dA =  km * As - (kp * (100-As) + kf * (100-As) * (100-As))
    dAs = kp * (100-As) + kf * (100-As) * (100-As) - km * As
    list(c(dA, dAs))
  })
}

times = seq(0,10,0.05)

out <- ode(y = yini, times = times, func = mich, parms = NULL)

plot(out) # this looks unstable

Michaelian system with ultrasensitive feedback

yini = c(A = 100, As = 0)
kp = 0.3; km = 1.5; kf = 0.1
kmf = 0.1; n = 3


mich = function(t, y, parms) {
  with(as.list(y), {
    dA =  km * As - ((kp + kf * (As^n/(As^n + kmf^n))) * (100-As))
    dAs = (kp + kf * (As^n/(As^n + kmf^n))) * (100-As) - km * As
    list(c(dA, dAs))
  })
}

times = seq(0,10,0.05)

out <- ode(y = yini, times = times, func = mich, parms = NULL)

plot(out)

Linear feedback plus saturating back reaction

kmb = 3
FR = function(As) {(kp + kf * As) * (100 - As)}
BR = function(As) {km * (As / (As + kmb)) * (100-As)}
plot(1:100,FR(1:100));lines(1:100,BR(1:100))

plot(BR(1:100))